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The properties of hydrogen-helium mixtures at Mbar pressures and 
intermediate temperatures (4000 to 10000 K) are calculated with 
first-principles molecular dynamics simulations. We determine the 
equation of state as a function of density, temperature, and compo- 
sition and, using thermodynamic integration, we estimate the Gibbs 
free energy of mixing, thereby determining the temperature, at a 
given pressure, when helium becomes insoluble in dense metallic hy- 
drogen. These results are directly relevant to models of the interior 
structure and evolution of Jovian planets. We find that the tem- 
peratures for the demixing of helium and hydrogen are sufficiently 
high to cross the planetary adiabat of Saturn at pressures around 5 
Mbar; helium is partially miscible throughout a significant portion of 
the interior of Saturn, and to a lesser extent in Jupiter. 

hydrogen-helium mixtures | high-pressure | helium solibility 

The two lightest elements, hydrogen and helium, are fas- 
cinating to physicists. Ubiquitous in the universe, their 
abundance ratio provide stringent checks on cosmological nu- 
cleosynthesis theories and the global distribution of hydrogen 
in the observable universe provides clues to the origin and 
large scale structures of galaxies. They are the essential ele- 
ments of stars and giant planets. Yet, in spite of the seem- 
ing simplicity of their electronic structure, there are many 
unanswered questions about their fundamental properties, es- 
pecially at high pressures. One such question is under what 
conditions are these elements miscible. The answer will have a 
crucial impact on our understanding of the evolution and the 
structure of the giant planets in our solar system and beyond. 

Jupiter and Saturn, the simplest among the Jovian plan- 
ets, are generally believed to have been formed approximately 
at the same time as the sun, although certain direct observa- 
tions (such as Saturn's excess luminosity) appear to contradict 
this planetary formation theory. In addition to being mostly 
made of hydrogen and helium, a characteristic of Jovian plan- 
ets is that they radiate more energy than they take in from 
the sun. Various models of their evolution and structure have 
been developed [1] |2j [3 ^ 4 to describe a relation between the 
age, volume, and mass of the planet and its luminosity. The 
current luminosity of Jupiter is well described with an evo- 
lution model for a convective homogeneous planet radiating 
energy left over from its formation 4.55 billion years ago. But 
a similar model seriously underestimates the current luminos- 
ity of Saturn [5]. Hence, either Saturn formed much later 
than Jupiter, or there is an additional energy source playing 
a more important role in Saturn than in Jupiter. In addition, 
the atmospheric abundance of helium in both Jupiter and Sat- 
urn appears to be lower than the accepted proto-solar values, 
more so in Saturn than in Jupiter [T]. 

Salpeter and Stevenson [6] [Jj El |9] proposed that helium 
condensation could be responsible for both the excess luminos- 
ity in Saturn and the helium depletion in the atmosphere of 
both Jovian planets. Suppose there is a region in the planet's 
interior where helium is insoluble; helium droplets will form 
and the denser helium will act as a source of energy, both 
through the release of latent heat, and by descending deeper 
into the center of the planet. Because Jupiter and Saturn 
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Fig. 1. Schematic phase diagram of hydrogen-helium mixtures at two composi- 
tions in the relevant range. Immiscibility lines (solid lines); this work: triangles (green 
- guide to the eye), Pfaffenzeller et al. 13 (brown), Hubbard-DeWitt 11 (turquoise), 
the latter two as parameterized in 0- In all three cases, the upper line corresponds 
to Xjj e = 0.0847 and the lower one to Xjj e = 0.0623. Isentropes (dashed 
lines); Jupiter assuming a composition Xn & = 0.07: squares (red - guide to the 
eye), Saturn assuming a composition X{j e = 0.0667: circles (blue - guide to the 
eye). Also shown in the lower left corner is the molecular H2 dissociation region of 
the mixture (dashed violet) from 19 and 20 . 



have different total masses, the thermodynamic conditions in 
the planetary interiors could be such that this condensation 
process is more prevalent in Saturn than in Jupiter. 

Previous attempts to calculate the immiscible tempera- 
ture, as a function of pressure and helium concentration, by 
Stevenson [7J, Straus [101, Hubbard, et al. [TT] Klepeis, et 
al. [12] and Pfaffenzeller, et al. |13| led to inconsistent con- 
clusions as to the importance of phase separation in the interi- 
ors of Saturn and Jupiter. The original theories of Stevenson, 
Hubbard and DeWitt were based on the assumption that the 
mixture consisted of fully pressure-ionized hydrogen and he- 
lium. For the temperatures and pressures found in Saturn 
and Jupiter, this assumption is now known to be inaccurate, 
especially for helium [141 115] . Klepeis, et al. [12] and Pfaf- 
fenzeller, et al. [13] developed mixture models based on den- 
sity functional theory (DFT). This opens up the possibility of 
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providing an accurate description of electron-ion interactions 
without assumptions on the extent of ionization. Klepeis et 
al. calculated the enthalpy of mixing at zero temperature 
from the analysis of crystal structures with different concen- 
trations of helium. Using those enthalpies and the assump- 
tion of ideal mixing for the entropy, they obtained a demixing 
temperature of 15000 K for XHe = 0.07, which suggests that 
there should be a major phase separation in both Jupiter and 
Saturn. However, this work neglected both the relaxation of 
the ionic crystal after the introduction of helium, and disor- 
der characteristic of a fluid. Using first-principles molecular 
dynamics (FPMD) simulations with the Car-Parrinello tech- 
nique, Pfaffenzeller et al. }13j developed a model including a 
realistic fluid structure. They performed molecular dynamics 
(MD) simulations of fluid pure hydrogen and estimated the 
free energies of a mixture by a re-weighting technique. They 
found a negligible temperature effect on the mixing free en- 
ergy up to temperatures of 3000 K and therefore disregarded 
thermal effects in enthalpies of mixing and used the ideal mix- 
ing for the entropy. They obtained immiscibility temperatures 
too low to allow for differentiation in either Jupiter or Saturn. 

In the present work, the temperature, pressure and com- 
position dependence of the enthalpy in hydrogen-helium mix- 
tures is computed with FPMD simulations (see Methods) 
based on DFT. DFT has become the method of choice in 
theoretical studies at high pressures, producing sufficiently 
accurate results for hydrogen and helium [161 117] . We neglect 
the zero point energy of the ions, which has been shown to 
be small and will have a negligible effect on the immiscibility 
temperature within our precision [18] . Using thermodynamic 
integration, we estimate the Helmholtz free energy of the mix- 
ture and determine the demixing temperature as a function 
of pressure and composition, thereby avoiding many of the 
previous assumptions and providing the most accurate pre- 
diction of the hydrogen- helium immiscibility, to date. Figure 
[I] summarizes the main findings of this work. The isentropes 
for Jupiter and Saturn determined from our DFT-based equa- 
tion of state (EOS) are shown along with the temperature of 
demixing. Overall, we find that the demixing temperature is 
high enough to support the scenario where helium is partially 
miscible over a significant fraction of the interior of the Jovian 
planets, with the corresponding region in Saturn being larger 
than in Jupiter. 



Results and Discussion 

We calculated the EOS of the hydrogen-helium system as 
a function of composition in the temperature range 4000 to 
10000 K and in the density range 0.3 to 2.7 g/cm 3 , by a series 
of FPMD simulations in the NVT ensemble. We studied 12 
different compositions to obtain an accurate interpolation of 
the energy and pressure. Using the EOS we calculated free 
energies by integrating along isotherms and isochores. We 
used the following multistep process to estimate the Gibbs 
free energy: 

1. We computed the free energy of an effective model at the 
"reference point" (T re/ = 10,000if, r s = 1.25) [] for all 
12 compositions (see Methods section). At this point the 
structure of the liquid mixture can be reasonably well re- 
produced by simple pair potentials between classical point 
particles. 

2. Using Coupling Constant integration (CCI) (see Methods 
section), we computed the free energy difference at the 
reference point between the DFT model and the effective 
model. 

3. Integrating the EOS along constant temperature and con- 
stant volume paths, we obtained the free energy difference 
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Fig. 2. Internal energy (o black), entropy contribution to the free energy (D red) 
and Helmholtz free energy (A blue) of mixing as a function of the helium number 
fraction, at the reference point. The red solid line represents the ideal entropy of 
mixing and the black and blue dashed lines are results from Pfaffenzeller et al. U3l . 



between the reference point and any other thermodynamic 
point in the range investigated. 
4. Finally, inverting the pressure-volume relations we ob- 
tained the Gibbs free energy of mixing as a function of 
pressure, temperature and composition. 

Figure [2] shows the Helmholtz free energy of mixing as 
a function of helium number fraction (xHe) at the reference 
point as well as its energetic and entropic contributions. The 
figure clearly shows, at least at the reference point, that the 
ideal mixing assumption to describe the entropy of mixing is 
very accurate for xh e < 0.2 but becomes less accurate for 
larger helium concentrations. The reason for this behavior is 
that the local environment of a proton in the low hydrogen 
concentration region is very different from the one it experi- 
ences in the metallic state of the pure system (see the discus- 
sion of pair correlations below). On the other hand, the inert 
character of helium makes it insensitive to change in the local 
environment in the low helium concentration region. In figure 
[2] we compare our results to the prediction of Pfaffenzeller et 
al. [13] who neglected thermal effects in the internal energy of 
mixing and used the ideal mixing law. The neglect of thermal 
effects in the internal energy results in a too large and negative 
mixing free energy. While it is true that the thermal effects 
are probably negligible at 3000 K, at this temperature the sys- 
tem is strongly immiscible so that the re- weighting procedure 
used in 43] to estimate those effects is likely to be inaccurate. 

In figure [3] we present the calculated Gibbs free energy of 
mixing as a function of composition; in panel a) several pres- 
sures are shown at a temperature of 8000 K, while in panel 



r s is the Wigncr-Scitz radius and defines the electronic density. 
2 We have concentrated the majority of our simulation efforts on the small 
XH e part of the phase diagram more relevant to planetary models. Quanti- 
tative prediction of miscibility at large xjj e is more difficult because of the 
smaller mixing free energy involved and will require additional investigations. 
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b) several temperatures are shown at a pressure of 10 Mbar. 
Note that at 8000 K, pressure has a small effect on the mix- 
ing free energy for low helium concentrations. In particular, a 
minimum in the free energy located at xhc ~ 0.1 is observed 
for all pressures investigated; this implies a stable mixture at 
this concentration. On the other hand, pressure has a strong 
effect at higher helium concentrations where, at a temperature 
of 8000 K, an increase from 4 to 10 Mbar eliminates the sec- 
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Fig. 3. The Gibbs free energy of mixing as a function of composition. Panel a), 
Data at 8000K for several pressures: 4 Mbar (□ - red), 8 Mbar (A - green), 10 
Mbar (V - blue) and 12 Mbar(o - violet). Panel b), Data at P= 10 Mbar for several 
temperatures: 5000K (□ -red), 7000K (A - green), 8000K ( V - blue) and 10000K 
(O - violet). Results of calculations are shown as symbols and error bars; solid lines 
are guides to the eye only. 
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Fig. 4. Demixing transition temperatures as a function of helium number frac- 
tion, for several pressures: 4 Mbar (□, black), 8 Mbar (A, red), 12 Mbar (o, blue). 
Continuous lines are the results of Ref. 1121 for 10.5 Mbar (brown) and Ref. 1131 for 
8 Mbar (green) and 10 Mbar (violet). 



ond minimum in the free energy curve|^] The common tangent 
construction is used to estimate the demixing temperatures. 
For points where no minima at high helium concentration is 
evident, we have assumed complete immiscibility. From the 
free energy plots of figure [3] it is clear that this assumption 
will have a negligible effect on the location of the minimum 
at low helium concentration. As shown in panel b) of figure 
[3j temperature has a strong effect on the mixing free energy, 
and hence, on immiscibility. An increase in temperature from 
7000K to 9000K (not shown in the figure) is enough to change 
the concentration of helium at the saturation point from 5% 
to 15%. 

Figure [4] shows the demixing temperature versus composi- 
tion for pressures ranging from 4 to 12 Mbar. Also shown are 
the results from the previous DFT-based calculations [121113] . 
As suggested by the free energy curve in figure [3] pressure 
has only a moderate effect on the immiscibility process. For 
a fixed helium fraction, the demixing temperature changes by 
approximately 500 K in a pressure range of 8 Mbar for the 
relevant concentrations (5% to 10%). 

Recent first-principles studies of pure helium have exam- 
ined the effect of temperature on band gap closure, suggesting 
that metallization in helium can occur at much smaller pres- 
sures than previously expected [H] . To examine the nature of 
helium in the mixtures, we calculated the electronic conduc- 
tivity of pure helium using the Kubo-Greenwood approach 
within DFTrl and obtained values well below 100 (ficm) -1 , 
even at the highest temperature and density reported here. 
Furthermore, in the recent work by Stixrude et al. [14], for 
p ~ 5.4 g/cm 3 the band gap is found to close at temperatures 
beyond 20000 K, well above our estimated demixing temper- 
ature. Metallization should enhance helium solubility, but as 
clearly shown here, for the pressures relevant to the modeling 
of Jovian planets, immiscibility occurs at temperatures well 
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Fig. 5. The effect of increasing concentrations of helium on the proton-proton 
radial distribution function at r s = 1.25 (top panel) and r 3 = 1.05 (bottom 
panel). 



These calculations were performed on 15 de-correlatcd configurations from 
a FPMD trajectory, where wc employed a 6x6x6 Monkhorst-Pack k-point 
grid. 
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below those required to produce ionization in helium [15] ; fully 
ionized models are not appropriate for describing the pressure 
dependence of the demixing temperature. At pressures much 
higher than those examined here, metallization of helium will 
play an important role and should produce significant changes 
to the pressure dependence of the immiscibility temperature. 

The structure of hydrogen is strongly influenced by the 
helium concentration. While at low Xn e hydrogen is in the 
mono-atomic fully ionized state, an effective proton-proton at- 
traction reminiscent of the molecular bonding develops upon 
increasing XHe, even at very high pressures and temperatures. 
Figure [5] shows several hydrogen- hydrogen radial distribution 
functions for mixtures with various helium concentrations, for 
temperatures of 8000 K and 10000 K and electronic densities 
given by r 3 — 1.05 and r 3 — 1.25 respectively. A molecular- 
like peak builds up smoothly as xh b —* 1- Under these con- 
ditions, helium is not ionized; this inhibits the delocalization 
of the hydrogenic electrons, enhancing the formation of weak 
molecular bonds. Because of the very low proton concen- 
tration, the observed proton-proton correlation can be inter- 
preted as resulting from an effective Morse potential. Fitting 
— \og(g pp (r))/T to this analytic form yields well depth pa- 
rameters ~300 times smaller than in an isolated hydrogen 
molecule. Such weak attraction gives proton pairs with short 
lifetimes, as also inferred from direct inspection of the MD 
trajectories. A similar stabilization of molecular hydrogen by 
helium, but at much lower temperature and density, has been 
previously reported close to the dissociation regime in pure 
hydrogen |19] . 

The computed demixing temperatures found here have im- 
portant implications for the study of the interior structure of 
hydrogen rich planets, especially Saturn. Our results sup- 
port the scenario where helium becomes partially miscible in 
the intermediate layers of the planet, with the excess helium 
falling towards the core through gravitational differentiation. 
This mechanism has been proposed to explain the high sur- 
face temperatures observed in Saturn [2] and the depletion 
of helium in the atmosphere. Whether the immiscible region 
will be large enough to explain all the observed properties of 
the planet remains an open question, but the current work 
represents a clear indication that a correct description should 
include phase separation. In general, DFT represents a large 
modification to the existing interior models, based mostly on 
the Saumon, Chabrier, Van Horn EOS [21] (SCVH). As re- 
cently shown by Militzer et al. [55], a combination of the DFT 
EOS with the inclusion of non-ideal effects in the mixing leads 
to isentropic Jupiter models that are cooler than the corre- 
sponding SCVH ones by about 1000-2000 K at high pressures. 
Even though the results reported here suggest that the inte- 
rior model can not be assumed to be isentropic, non-adiabatic 
models based on DFT including non-ideal effects should pro- 
duce qualitatively similar results to those found by Militzer. 
In order to give an indication of the general changes expected 
with the proper inclusion of these effects, we calculated the 
SCVH isentropes (corresponding to the interpolated EOS) of 
Jupiter and Saturn using our DFT-based EOS. The results 
are shown in figure [l] where we have assumed fixed helium 
concentrations of XHe = 0.07 for Jupiter and XHe = 0.0667 
for Saturn. 

At low density (p < 0.3 g/cm 3 ), before the dissociation of 
hydrogen plays a significant role, the pressures produced by 
the SCVH model are in good agreement with DFT. In order to 
compare the entropy of both models at low density, and assess 
the agreement of the isentropes there, we performed simula- 
tions with XHe = 0.0667 following a path through the region 
of miscibility starting from the reference point and ending at a 
density of p = 0.3 g/crn 3 , where we simulated several temper- 

4 I 



atures. Then using thermodynamic integration, we calculated 
the free energy and the entropy at the point where the SCVH 
planetary isentrope crosses this density. We obtain entropies 
that agree with the SCVH model to within 1.5%. With the 
onset of dissociation for higher densities, the isentrope from 
the SCVH model deviates significantly from DFT results. Al- 
though the possible immiscibility puts in question the use of 
an isentropic model for the interior of Saturn, our analysis pro- 
vides an estimate of the magnitude of changes to the SCVH 
based models. The present results are in good agreement with 
those reported by Militzer 22 for the Jupiter isentrope. As 
can be seen, the portion of the interior of Saturn correspond- 
ing to pressures between 1.5 to 5.5 Mbar should be inside the 
immiscible region, with the stable concentration of helium de- 
pending on pressure. 

Conclusion 

In summary, we have carried out an extensive investigation of 
the properties of hydrogen-helium mixtures at pressures and 
temperatures that are relevant to the interior of Jupiter and 
Saturn using state-of-the-art ab initio simulation methods. 
By using a combination of first-principles molecular dynam- 
ics simulations within DFT and thermodynamic integration 
techniques, we have accurately determined the Gibbs free en- 
ergy of mixing over a wide range of density, temperature and 
composition. Our work differs from previous investigations in 
that it does not rely on any assumptions about mixing func- 
tions. Our simulation results are consistent with the idea that 
a large portion of the interior of Saturn has conditions such 
that hydrogen and helium phase separate; this can account 
for the apparent discrepancy between the current evolution 
models for Saturn and observational data. 

The accuracy of the present results is primarily limited 
by the approximate density functional used. Quantum Monte 
Carlo calculations for helium and hydrogen do not require 
assumptions about the electron correlation, pseudopotentials 
and zero-point energy of the nuclei. A systematic investiga- 
tion of the correction to DFT using the Quantum Monte Carlo 
method described in |17| is in progress. 

Methods 

First Principles Calculations. The FPMD simulations per- 
formed in this work were based on Kohn-Sham den- 
sity functional theory using the Perdew-Burke-Ernzerhof 
(PBE) exchange-correlation functional. We employed Born- 
Oppenheimer MD (BOMD) within the NVT-ensemble (with 
a weakly coupled Berendsen thermostat), as implemented in 
the Qbox code (http:/ /eslab. ucdavis.edu/software/qbox). We 
used a Hamman type [23] local pseudopotential with a core 
radius of r c — 0.3 au to represent hydrogen and a Troullier- 
Martins type |24j nonlocal pseudopotential with s and p chan- 
nels and r c — 1.091 au to represent helium. Tests have estab- 
lished the accuracy of the hydrogen and helium pseudopoten- 
tials over the relevant pressure range, up to 13 Mbars. ^ A 
plane wave energy cutoff of 90 Ry was used for r s > 1.10 and 
of 115 Ry for r 3 < 1.10. Empty states were included with an 
electronic temperature set to the ionic temperature. To inte- 
grate the equation of motion during the dynamics we used a 
time step of 8 a.u. 



4 The sensitivity of the free energy to the choice of the helium pseudopoten- 
tial was estimated by comparing the enthalpy of 20 ionic configurations with 
a Hamman type local pseudopotential with r c — 0.218 au for helium and a 
plane wave cutoff of 450 Ryd. The difference in the free energy of mixing is 
estimated to be 0.1 mH/atom, this is smaller than the error bars reported in 
Fig. [3] 
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We used 250 electrons for r s < 1.65 and 128 electrons 
otherwise. The Brillouin zone was sampled at the F-point. 
In order to reduce systematic effects and to get accurate pres- 
sures, we added a correction to the EOS designed to correct for 
the plane wave cutoff and the sampling of the Brillouin zone. 
To compute this, we studied 15-20 configurations at each den- 
sity and composition by using a 4x4x4 grid of k-points with a 
plane-wave cutoff of at least 300 Ry. The actual plane-wave 
cutoff used depended on density and was chosen to achieve 
full convergence in the energy and pressure. 

For the calculation of the EOS, we studied 4 temperatures, 
6 densities and 12 compositions for a total of 288 simulations. 
We also performed 15 simulations to extend the calculation of 
the free energy to low density. For the integration of the free 
energy we studied 8 compositions, each one required 5 addi- 
tional calculations. In the case of the BOMD simulations, we 
first equilibrated the system using a suitable effective model 
and subsequently allowed 300-500 timesteps of equilibration 
with DFT, averages were accumulated for approximately 2000 
time steps. For the thermodynamic integration, the simula- 
tions were first equilibrated with the effective potential and 
subsequently allowed to equilibrate for 1000-2000 timesteps 
using the mixed DFT-effective potential, averages were cal- 
culated for 12000-18000 timesteps. In both cases, this was 
sufficient for accurate results. In total, the simulations re- 
ported here used approximately three million CPU hours on 
a large Opteron-based Linux cluster. 

Coupling Constant Integration. CCI allows us to calculate the 
difference in free energy between systems with different inter- 
acting potentials. For a system described by the potentials Vi 
and V%: 



V(X) 

F 1 (T,V,N)-F 2 {T,V,N) 



XVi + (1 - X)V 2 
1 dF(X) _ 

o 



'-dX 



(V, - v 2 ), 



[1] 



d\, [2] 



where {) TVNX represents a canonical average with the po- 
tential V(X). Any functional form of the two potentials is 
formally allowed in Eq. [I] but the use of similar potentials 
makes the integration of the free energy difference consider- 
ably easier in practice. 

To represent the interaction between the atoms in the clas- 
sical system, we used reflected Yukawa pair potentials: 



V{r) = 



+ 



-HL-r) 
(L-r) 



.-bL/2 , 



r < L/2, 
r > L/2, 



where a,b and L are free parameters and depend on the iden- 
tity of the atom^] As shown in figure [HJ this potential was 
found to exhibit similar pair correlations as the DFT model for 
pure hydrogen, but was not as good for helium. We choose the 
pair potentials such that the effective model was fully misci- 
ble at the reference point to avoid crossing a phase line during 
the integration. We computed the Helmholtz free energy of 
the effective model using CCI and classical MC simulations, 
with the second potential set to zero in Eq. [T] From this 
we determined the free energy of the effective model since the 
free energy of the non-interacting model is known. We also 
calculated the free energy by integrating the pressure from 
the reference volume to a volume large enough that the sys- 
tem is ideal; the pressures were obtained by a series of classical 
MC simulations. Both approaches produced agreement within 
noise. 



The free energy difference between the DFT-based and the 
effective models was calculated using the CCI approach. Fig- 
ure [6] shows a comparison of the radial distribution functions 
of the effective and DFT models for selected compositions at 
the reference point. To compute the required canonical aver- 
ages, we used a Hybrid Monte Carlo (HMC) algorithm |25] . 
which allows for an efficient sampling of large systems with 
many-particle moves. HMC results in exact sampling of the 
canonical ensemble without time step errors. In this work, the 
HMC approach was found to be as efficient as MD if the time 
step is chosen carefully. Figure[7|shows the results of the HMC 
simulations for several compositions. The curves are smooth. 
This is the only requirement to justify the procedure. Fig- 
ure]^ summarizes the main results of the computations at the 
reference point. 
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Fig. 6. Compa rison of several radial distribution functions between the effective 
(blue) and the DFT (red) system at the reference point. We choose Yukawa pair 
potentials to calculate the Helmholtz free energy at the reference point. 
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7. Results of the simulations for the mixed effective-DFT potential used 
in the CCI of the free energy difference for several compositions. As can be seen, 
the results depend smoothly on the coupling parameter and no singular behavior is 
observed close to the endpoints. 
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